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Abstract. We describe and analyze an interior-point method to decide feasi- 
bility problems of second-order conic systems. A main feature of our algorithm 
is that arithmetic operations are performed with finite precision. Bounds for 
both the number of arithmetic operations and the finest precision required are 
exhibited. 



1 Introduction 

It is now widely accepted that the most efficient algorithms for solving the gen- 
eral type of second-order conic problems are interior-point methods (IPMs). IPMs 
infallibly demonstrate very fast numerical convergence, by far outperforming their 
theoretical estimates. 
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Second-order conic programming problems contain linear programming problems 
as a special case, and at the same time can be embedded into the class of semidefinite 
programming problems. It is, however, not advisable to solve SOCP problems by 
semidefinite programming methods (see [1], [15]) as IPMs that solve SOCP directly 
have a much better complexity (both in theory and in practice). SOCP problems 
have lately received considerable attention due to their many applications [15]; they 
appear to be at the boundary of the problems for which interior-point methods 
can solve large instances, a fact that is linked to the implementation of commercial 
software for the solution of second-order programs such as MOSEK 1 or CPLEX 2 . 

We are interested in solving homogeneous second-order conic feasibility prob- 
lems. That is, given a second-order cone K C H n and a matrix A 6 IR mxn , decide 
which one of the primal-dual pair of problems 

Ax = 0, A T y + s = 0, 

x hx 0, 1 ' s^kO, { ' 

is strictly feasible (i.e., the relevant conic constraint is strict) and provide a solution 
to the feasible problem. It is well-known that each of (P) and (D) above has a strict 
solution if and only if the other one has no nonzero solutions. 

Recall that a second-order cone is a direct product of a finite number of Lorentz 
cones. The Lorentz cone C p C R p+1 is defined to be 

C p = {x G R p+1 | x > ||x||}, 

where for a vector x G R p+1 indexed from to p we let x = (x\, x%, . . . , x p ) G W. 
For our primal-dual pair of problems (P)-(D) we take 

K = C>ni X C n2 X • • • X C nrl 

where r is the number of the Lorentz cones comprising K, and J2l=i n % = n with rij 
being positive integers for all i from 1 to r. 

We propose a finite-precision algorithm for solving the SOCP feasibility problem 
and provide rigorous bounds for the finest machine precision and the maximal num- 
ber of iterations needed. The proposed algorithm is designed to work with variable 
precision, that is, the machine precision can be re-adjusted along the way. 

Our bounds depend on Renegar's condition number [12], [13], which is consistent 
with similar bounds obtained for the polyhedral case in [5]. Let pp(A) and Pd(A) 
be the distance to infeasibility of (P) and (D) respectively defined by 

p P (A) = inf{||AA|| : (A + AA)x = 0,x >- K is infeasible} 

and 

p D {A) = inf{||AA|| : -(A + AA) T y >~k 0,t/ £ IR m is infeasible}. 

1 http: //www.mosek. com/ 

2 http: //www. ilog. com/products/cplex/ 
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Renegar's condition number C (A) is defined as the reciprocal of the relative distance 
to ill-posedness of the pair (P)-(D): 



max{pp(A),/9 D (A)}' 

Although any equivalent matrix norm can be used to define C {A) , in our analysis we 
choose to use the standard operator norm induced by the Euclidean scalar product. 
We say that the problem is ill-posed if both pp(A) = pd(A) = and hence C(A) = 
oo. 

Our main result, Theorem 1, shows that there exists a finite precision interior- 
point method which, with input a matrix A £ H ra and a second-order conic structure 
K (consisting of r Lorentz cones), decides which one of the two systems (P) or (D) 
is feasible. We estimate both the number of iterations of the algorithm and the 
precision required as functions of the size of the matrix, the number r of Lorentz 
cones in K, and the condition number of the problem. The finest required precision 
is 

1 

II = 

0{{m + nf/ 2 r l ^C{A)V^y 

and the number of main interior-point iterations performed by the algorithm is 
bounded by 

0(r 1/2 (logr + \ogC{A)). 

Strictly speaking, our algorithm solves both the decision problem — decide which 
one of the problems (P) and (D) is feasible — and the function problem — if either 
one of the problems (P) or (D) is strictly feasible produce a (possibly approximate) 
solution for it. 

Throughout the paper, we use standard notation wherever possible. We in- 
dex our variables according to the second-order conic structure. That is, x = 
(x±, X2, ■ ■ ■ , x r ), s = (si, S2, ■ ■ ■ , s r ), where Xi,Si G R ni+1 for all i = l,...,r. 
Throughout the paper we assume that ||^4i||f = l/v^> an d hence ||^4||f = 1- Note 
that this assumption is trivial from a computational viewpoint; if Ai / mxni , it 
takes a few operations to reduce the matrix to this form and it is easy to recover 
solutions of the original system from those for the reduced one. The condition num- 
ber of the new matrix may change, however. But one can show as in [5, §11.2] that 
this change can not be large. 

Finite precision analyses are pervasive in Numerical Linear Algebra; they are 
much less common in optimization. While the effects of finite precision when solv- 
ing linear programming problems had been early noticed (e.g. [2, 4, 8, 14, 16, 20]) 
there was no condition-based round-off analysis even for linear programming prob- 
lems until recently. This was done for the feasibility problem for polyhedral conic 
systems [5], for the optimal value of linear programs [18], and for the computation 
of optimal basis and optimal solutions of linear programs [3]. To the best of our 
knowledge, our work is the first such analysis for nonlinear cones. 
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Our paper is organized as follows. In Section 2 we use a relaxation scheme 
introduced by Peha and Renegar [10] and Vera et al. [19] to reformulate the feasi- 
bility problem via an optimization one and remind the basic idea of interior-point 
methods. Then we relax the standard results of IPM analysis to make room for 
computational errors. We do not deal with finite-precision issues directly until Sec- 
tion 3, where we describe our algorithm in detail and estimate errors arising on every 
step of floating-point computations. The last section is devoted to the proof of the 
main result, and essentially fits the error estimates obtained in Section 3 into the 
gaps made to this purpose in our extension of the IPM analysis done in Section 2. 



2 Interior-point method for SOCP feasibility problem 



We use a relaxation scheme introduced by Peha and Renegar in [10] and later 
extended in [19]. This relaxation scheme reformulates the feasibility problem (P)- 
(D), for the more general case when K is a symmetric cone, as a pair of primal- 
dual optimization problems in higher dimension and solves this pair by a standard 
short-step interior-point method. We next summarize the main ingredients of this 
approach. 

It was shown in [19] that the pair (P)-(D) is equivalent to the following primal- 
dual pair of optimization problems 



■ —T - 

mm c x 



s.t. Ax = b 
x^kO 



(P') 



and 



— *-p 

max b y 
s.t. A y + s = c 
shK 0, 



(D') 



where K, = K x C n x C m , x = (x, t, x', r, x"), s = (s,t s , s' ,r s , s") £ 
y = (y,y',-r ] )€M m+n+1 , 
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This equivalence should be understood in the following sense: If p(A) > then 
a primal-dual interior-point method applied to the pair (P')-(D') yields a strict 
solution to whichever of (P) or (D) is strictly feasible. In particular, since the 
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optimal value of the pair (P')-(D') is zero, a corresponding strict solution to the 
original feasibility problem can be straightforwardly recovered from the first entries 
of the extended variables (x in the case of (P) and y and s in the case of (D)). 

In the sequel, to simplify notation, we will denote m := m + n + 1 and n := 
2n + m + 2 so that A G R mxn , b G IR m and c G R n . We also let r := r + 2 and 
V = {{x,y,s) | x G int(/C), y € R m ,se int(/C)}. 

The pair (P')-(D') can be solved via a primal-dual interior-point algorithm. We 
refer the reader to [13] and the references therein for a detailed exposition of the 
theory of IPMs. We next recall the concepts and results from this basic theory that 
will be used in the paper. Consider the following self-scaled barrier function for the 
cone K.: 

= ~ (E ln 04 " \\^f) + Ht 2 - ||x|| 2 ) + ln(r 2 - ||z"|| 2 )^ . (2.1) 

Let g(x) = Vf(x) and H{x) = V 2 f(x) denote respectively the gradient and the 
Hessian of /. Let e G K, denote the unique point such that H(e) = I, that is, 
e = (e±, . . . , e r ) where e^o = 1 and = for i = 1, . . . , r. 

Sometimes we will also need to work with the self-scaled barrier function for the 
cone K: 



/(z) = -J>( a 

i=l 

and will let g(x) = V/(x) and H(x) = V 2 J(x). Notice that H{x) 



H{x) J) 
H{x)Y 

where H(x) denotes the Hessian of the function — ln(i 2 — ||x'|| 2 ) — ln(r 2 — ||x"|| 2 ) 
Given x G /C, the local norm \\ ■ \\% in IR m is defined as 



:= u T H(x)u. 

Likewise for x G K. 

The central path of (P')-(D') is the set of solutions {(x(fi),y(fi), G V : n > 
0} to the system of equations 

Ax = b 

A T y + s = c (2.2) 
s + fig(x) = 0. 



Given z = (x, y,s) G V define 



X L S 



2r 

Note that if z belongs to the central path for a certain value of ^ then fi(z) = [i. 
We may sometimes write fi for ^(z) when z is clear from the context. 
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The basic idea of a path- following interior-point method is to generate a sequence 
of points on a suitable neighborhood of the central path that converges to optimality. 
The suitable neighborhood is the following. 



Definition 1 Given (3 G (0, 1/15), the central neighborhood Np is defined as the set 
of points z = (x, y, s) G V, such that the following constraints hold: 



Ax = b 
A T y + s = c 

S + V(z)g(x)\\-U(z)g(x) < P- 



The main computational step of each interior-point iteration is to solve a lin- 
earization of the central path equations (2.2) at the current iterate z 6 Mr. The 
linearization that we will rely on is as follows. We will see (cf. Proposition 1(f) 
below) that for all x, s G int(/C) there exists a unique scaling point w G K. such that 



Given z = (x,y,s) G Np, the Nesterov-Todd direction (Ax, Ay, As) is the solution 
to the following linearization of (2.2): 



By [19, Proposition 4.6], the initial point in step (i) of Algorithm IP below is 
in the central neighborhood Mp. Furthermore, by [19, Propositions 4.4 and 4.5] if 
the original pair (P)-(D) is well-posed (i.e. C(A) < oo), then a point z G Mp with 
fi(z) small enough yields a strict solution to either (P) or (D), whichever is feasible. 
Indeed, by [19, Theorem 3.1] Algorithm IP halts in at most 0(y / r(logr + logC( J 4))) 
iterations and yields a solution to either (P) or (D). (See [19] for details.) 

Algorithm IP (A) 

Let f3, 5 G (0, 1/2) be fixed constants such that 



H(w)x = s. 



AAx 
A T Ay + As 
Ax + H(w)~ 1 As 






-(fj-g(s) + x). 



(2.3) 




(i) Let 



a := — ; M : = 



a||^4e 



and 



x = (cue, 1, ae, 2M, —aAe) 



* = (^,-M e ,l,0). 
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(ii) If A T y < K 0. then HALT and 

return y as a strictly feasible solution for (D). 

(iv) If a miQ (AH(x)~ 1/2 A T ) > rfi(z), then HALT and 

return x + H(x)- 1 A T (AH(x)- 1 A T )- 1 strictly feasible so- 

lution for (P). 

(v) Set 77:= (l-^)^). 

(vi) Compute Az := (Ax, Ay, As) by solving (2.3) for \x = fx and up- 
date z by setting 

z + := z + Az 

(vii) Go to (ii). 

It should be noted that the analysis in [19] assumes that all computations are 
performed with infinite precision. Our initial step for a finite-precision algorithm is 
to show that the results in [19] can be extended to make room for computational 
errors. In particular, Lemma 1 below shows that even if the system (2.3) is solved 
inexactly, we can still ensure that the iterates remain in the central neighborhood. 

Lemma 1 Let z G Mp, ~p = (1 — ^)fi(z) with \5' — S\ < ^ and z + = z + Az be 
such that 

AAx =0 

A T Ay + As =0 (2.4) 
Ax + H(w)~ 1 As = -(Jigis) + x) + r, 

Where INI ^ 120r(2 A 'r ( / !(l)+l) ' Then Z+ G M P and \^ Z+ ) ~ Z 1 ! ^ l20F^ 

PROOF. See Section 2.2. □ 

We note that when the solution Az to (2.3) is computed exactly, the point 
z + := z + Az satisfies ^(z + ) = ~p. For details, see [13]. 

The following two lemmas are in the same spirit as [19, Propositions 4.4 and 
4.5]. In particular, they guarantee that if either (P) or (D) is strictly feasible then 
a point z G Np with u.(z) small enough yields a strict solution to either (P) or (D). 
Lemma 2 provides the relevant bound for /j,(z) in the case (P) is strictly feasible and 
Lemma 3 does so for a strictly feasible (D). 

Lemma 2 Let z = (x, y, s) G Np and assume pp(A) > 0. Then 

a min (Aff (*)-^ T ) > ( {1 ~flZ iA) ) 2 " W) 2 - ( 2 - 5 ) 
The latter in turn implies that if fi(z) < ^f^+pf^ then the point 

x + H(x)~ 1 A T (AH(x)- 1 A T )- 1 x" 
is a strict solution to (P). 
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Proof. This is an immediate consequence of the proof of Proposition 4.4 in [19, 
pages 259-260]. □ 



Lemma 3 Let z = (x, y, s) G Np and assume Pd(A) > 0. Then 

\\H(s)\\ < 
In particular, for i = 1, . . . , r 



4r 2 



sm ~ \\si\\ > ^-j=p D (A). (2.6) 

Furthermore, if p(z) < 4r 2 (l — /3)p£>(A) then y is a a strict solution to (D). 

Proof. This is an immediate consequence of the proof of Proposition 4.5 in [19, 
page 260]. □ 

The rest of this section is devoted to proving Lemma 1 and a technical lemma 
related to the conditioning of the matrix arising at each interior-point iteration of 
Algorithm IP. In §2.1 we state and prove a few technical statements which will be 
employed in the subsequent proofs. Then in §2.2 we prove Lemma 1. The proof 
of Lemma 1 is a straightforward adaptation of the proof of Theorem 3.7.3 in [13]. 
Section 2.3 presents Lemma 11, which is similar in spirit to Lemma 2. This technical 
result will be crucial in our finite precision analysis in Section 3. 



2.1 A few useful relations 

The analysis of IPMs heavily relies on the properties of the barrier function. Here 
we briefly remind a few essentials that will be used later. More details can be found 
in [13]. The barrier function / gives rise, for each point x in the domain Df of /, 
to a local inner product ( , ) x induced by x and defined by 

(u,v) x = (u,H(x)v). 

The local norm || || x is then given by \\v\\ x = (v,v) x 2 . In the local inner product 
( , ) x , the gradient at y is g x {y) '■= H(x)^ 1 g(y) and the Hessian is H x (y) := 
Hix^Hiy). 

Our function / defined by (2.1) is a self-scaled barrier with the barrier parameter 
v = 2r. We will also use single components of /: f% = — h^x 2 ^ — ||xj|| 2 ) for all 
i = 1, . . . , r. For each /j the barrier parameter is v = 2. Our development relies on 
the following key properties of self-scaled barrier functions [13]. 

Proposition 1 Let f be a v ■-self-scaled barrier function and x G Df. 
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(a) If \\y — x\\ x < 1 then y G Df and, for all v 7^ 0, 

1 - \\y - x\\ x < 1% < ! 



1 - lb - z|U' 

(b) {z£D/:(z - x,g(x)) > 0} C {z : ||z - ar|| x < z/}. 

(c) -g(x) G = x, H{-g{x)) = H(x)~\ and \\H{x)- l \\ < \\x\\ 2 . 

(d) Fora>0 

g(ax) = —g(x), and H(ax) = —^H(x). 
a cr 

(e) H{x)x = — g(x) and (x,g(x)) = —v. 



(f) Given another point s G there exist a unique "scaling point" w G Df such 
that 

H(w)x = s, and H(w)g(s) = g{x) 
and a unique "reverse scaling point" w* G Df such that 

H(w*)s = x, and H(w*)g{x) = g(s). 

Furthermore, w* := —g(w) and, for all (i > 0, the points w := ^/Jlw and 
w* := sfjlw* satisfy 

\\x 



and 

,1 4 

\s + ng{x)\\^ x ) > mm 



\\y-*\\ 



(g) If \\y - x\\x < 1 then \\g x {y) - g x (x) - H x {x){y - x)\\ x < T 

(h) If\\x-y\\ y <l then\\v\\_ gy(x) <(l + \\x-y\\ y )\\v\\y for allv £TR n . □ 

Lemma 4 Let z G J\fp and denote ji = fi(z). Then \\x\\ = \\x'\\ < 1, ||x"|| < r < 2r/i, 
< 1, \\s'\\ <t s < 2r/x, and \\s\\ < 2r/z + 1. 

Proof. The bounds on x, x' and s" follow from the equalities Ax = b and 
A^y + s = c together with x, s G /C. The inequalities ||a;"|| < r < 2r/i and ||s'|| < 
t s = rj < 2v\i follow from the equalities r + rj = c^x — b^y = x T s = 2r/x. Since 
z G N/3, we have s = As" - s' and therefore ||s|| < ||s'|| + ||A||||s"|| < 2r/i + l. □ 

The next lemma bounds the norm of the scaling matrix using the results above. 
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Lemma 5 Assume z € Np, with (3 < 1/15. Then 

ii^w-n< 4+9 '; ( f r2 < (2+3 f ( : )r)2 

H(z) n(z) 

and 

PROOF. By Proposition 1(f) we have — w\\w < f/3, therefore, by Proposi- 
tion 1(a) and (c), respectively, 



II K ) II- (1 _ 5^2" ^ (l-f/3)2 H II 

From Lemma 4 we have ||x|| 2 < 3 + 8u.r 2 (with /i = fi(z)), hence, 

II rrf \~H X \\tt(-\-H <r 3 + 8 A 2 , 4 + 9A 2 

II = -I^H 11^ (T^l^jv^ — ^ — • 

Similarly, applying Proposition 1(f) and Lemma 4 to H(w*)~ l we have 



W)-'n<- l|sT ' 4(2rf,+1)2 

A* ( 
Lemma 6 Let z e J\fp with (3 < ^ . Then 



iifhii = -h^ti < n 11 m < v 7 y • □ 



u(z) 1/2 

\\H(w)-^(-pg(x) + m<^^- 

Proof. Since z G JV^, we have from Proposition 1(f) that 

4 4 

P> + s||- w (f) = gp-^l|w= - sib*- ( 2 - 7 ) 

Then, applying Proposition 1 and (2.7) two times and using (3 < 1/5, we obtain 

llw/M + ^ 77-5^ — < (1 _ 5p)( 1 _p) ^ pMfr) + s||- w (f)- (2-8) 

By the triangle inequality and (2.7) 

+ *ll-^) < ^-^\\g(x)\\- 9 (x) + + ^ll- W (a) < ^- (2-9) 

From (2.7) and (2.9) we have 

\\H(w)-^(-pg(x) + s)\\ = » 1/2 \\Mvx) + ^ < □ 
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Lemma 7 Let z G Np for some P<\. Then fori = I,..., r 

xJ Si > 2(1 - p)n (2.10) 

and 



PROOF. From ||sj + W^OII-^g^ ) — <^ 2 an< ^ Proposition 1 we have 

/3 2 >IN| 2 _^)--<^} + 2. (2.12) 



(4 " IN! 2 )(4 " INI 2 ) > 4(1 - /3)V 2 - (2-11) 

2 

From Proposition l(a,e) 

llsill-^0 > INU (l " II* + M(xi)\\-M( Xi) ) > V2(l - (3). (2.13) 
Now (2.10) follows from (2.12) and (2.13). By Proposition 1 

< \\si + ng(xi)\\ 2 s . = \\si\\ 2 s . - 2fi(g(x),g(s)) + || - fig(xi)\\ 2 s . 
and by the definition of / 

l t \ i \\ ^\ x ii s i) 

ig{xU{s)) = (4-MP)(*-MI')- 

Therefore, 

{xl " IN|2)(4 " IN|2) - w^uxiWs; (2,14) 

By Proposition 1 

n-^,)ii„< "-'f-; 11 -^-' <-4. ( 2 .i 5 ) 

Now (2.11) follows from (2.10), (2.14) and (2.15). □ 
Lemma 8 Let z G Np. Then 

t>(1-/3)u.(z). (2.16) 

Proof. By Lemma 7, taking i = r + 1 in (2.10), and using ||s"|| < 1 from 
Lemma 4 

2(1 - P)n < tt s + x" T s" < r(l + ||s" ||) < 2r, 
which yields (2.16). □ 
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2.2 Proof of Lemma 1 

Let w be a scaling point of the pair x, s. Then by Lemma 5 and the assumption on 
the norm of r we have the following bound 

IHI- < J-||tf M 1/2 IIIIHI < J- 2(2r/i + 1) £ <±- (2 17) 

\\r\\ w S ^ 1/2 \\M[w) ||||r|| < ^ ^ 1/2 i2 (2r/i + l) " 60' 10 
Since z G A/"/?, and 5' < ^ • §§ < ^§ < , we have 

\\s + ~Pg(x)\\g{x) < \\S + fJ,g(x)\\ g (£) + |//-/Z| \\g(x)\\- g (£) 

Therefore, by Proposition 1(f) 

5 5 15 

\\x-w\\w< j\\s+~pg(.x)\\_ Jlg{£) < - ■ — = — . (2.18) 

Define 

u := gw{x) + 2W - x = gw(x) - gm(w) - H w (w)(x - w). ( 2 -19) 
Since \\x — w\\yj < ^ < 1, Proposition 1(g) yields 

\X — . T7T)4 25 

52 

Observe that 



u\\ w = "", <^Z04 = ^_. (2.20) 

1- i-iL ~ 1- A 52-47 



Ax + H(w)~ 1 As = -H(w)~ 1 (s + -pg(x)) + r 

= —x — JiH(w)~ 1 g(x) + r by Prop. 1(f) 

= —x — H(w)g(x) + r by Prop. 1(d) 

= -x- g w {x) + r 

= 2(w-x)-u + r (by (2.19)). 

Hence w — x = ^(Ax + H(w)~ 1 As + u — r) and so 

w — x + = ^(—Ax + H(w)~ 1 As + u — r), and 
w - H(u>y 1 s + = ^(Ax - H(w)~ l As + u -r). 
Consequently, 

H(w)~ l s + = 2w — x + — u + r. 
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Since Ax _L„, H(w) As, we have 

< 2 (II ~~ Ax + H(w)~ 1 As\\ w + \\u\\ w + \\r\\ w ) 

= ^(\\Ax + H(w)- 1 As\\ w + \\u\\ w + \\r\\ w ) 

._ _ 1 1 „ 1„ , 1„ 
= \\w - x - -u + -r\\ w + -\\u\\ w + -\\r\\ w 

— 11^ ^HuT - 1 - 1 1 1 1 u7 ~ l~~ ll^lli/J 

^ 25 1 fi 

" H2 + 52^7 + 47 = 47 to P-17), (2.18), and (2.20)). 



Thus 



||fr(w) L s + +g w (x + )\\ w = \\g w {x + ) + 2w - x+ - u + r\\ w 

= Wgwix^) + 2W - x + \\ w + \\u\\w + \\r\\w 

= \\9w{x + ) ~ 9w{w) - H w {wY l {x + -w)\\ w + \\u\\ w +\\r\\ 



I — -M-II2 
\W — X^ U " 



< z rr= ^TTT-II + \\ u \\w + \\r\\w (by Prop. 1(g)) 

1 — \\W — X + \\w 

36 25 1 107 . „, 

£ 47-41 + 2444 + «i < 2132 (M2.17) and (2.20)). 



From Proposition 1(h) we have for all v € M n 

(6 \ 53 
1 + M\w = ZZz\\v\\w- 

Thus 

\\s+ +~fig{x + )\\-g^+) = -fl\\H{w)- l s^ + gw{x + )\\- 9w (5;+) , oU 

< 53 . ML- a< (± -001) -n ^- ZL > 

v 47 2132^^ 1 15 /i. 

To finish we need to show that fJ-(z + ) is close to /Z. By our assumption we have 

Ax + #(i/;) _1 As + x = —Jlg(s) + r. 
Taking inner product with s and using Proposition l(e,f) we get 

(s, Ax) + (x, As) + (s, x) = -7l(s, g(s)} + (s, r) = 2r7i + (s, r). (2.22) 
Since Ax _L As, we have (x + Ax, s + As) = (x, s) + (x, As) + (Ax, s), so by (2.22) 

M ( z +) = J_(f + Ax, s + As") = ((x, s) + (x, As) + (Ax, s)) = 71 + (2.23) 

Using (2.23), Lemma 4 and the assumption on r, we get 

_7i| - Kfl^ < _J^_ - E < A. ( 2 24) 

} " 2r £ 12Qr> " 12Qr2 (l - < 120" ^ 
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Therefore, from (2.21) and (2.24) we get 

+ ^z+)g(x+)\\_ gia;+) < 11^+^(^)11-^+) + \^z + )-n\w + )\\4m) 

= Us 1 " + ~P-g{x r )\\- g (x+) + 2r\fi(z + ) -pi 

£ (i- - 01 ) F+ ^ < ^" (2+) - 

We have already shown that ||W — x + ||w < J=, therefore by Proposition 1(a) 
x + G /C. Since by Proposition 1(c) this yields —fi(z + )g(x) G /C, from 

and by the same argument we have G /C, and hence 

z+ G £>. (2.26) 

Observe that z + satisfies the linear equations Ax + = 0, X r y f +s f = by definition. 
Together with (2.25) and (2.26) this yields z + G Hp. □ 

2.3 On the condition of the matrix AH(wy l/2 A T 

The purpose of this section is to present Lemma 11 which will be crucial in our 
finite precision analysis in Section 3. This lemma is in the same spirit as Lemma 2. 
While it will not be used until Section 3, we choose to place it here to facilitate 
understanding of the proof. 

We will rely on the following key characterization of the distance to ill-posedness 
due to Renegar [12, Theorem 3.5]. 

Proposition 2 For any given linear operator A : R n — > R m and any cone K 
p P {A) = sup{<5 : ||t>|| < 5 => v G {Ax : ||x|| < l,x G K}} 

and 

p D (A) = sup{5 : \\u\\ <5^ue {A T y : \\y\\ < 1} + K*}. □ 

We will also rely on the following perturbation result, an extension of [19, Theorem 
5.1]. 

Lemma 9 Let f3 < yg, z G Afp. Assume b G lR m is such that (y,b) < and 

1 



Av = b 



> 1. 
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If a > then the optimal value of the perturbed problem 

min (Fu 



s.t. Au = b + ab, (2.27) 



is greater than (Fx. 



PROOF. This readily follows by putting together [19, Theorem 5.1] and Proposi- 
tion l(a,f). □ 

Lemma 10 Let Ab = (Aft 7 , Ab 11 , Ab 111 ) G M m+n+1 = R m . Then there exists u 
such that 

A T u = b + Ab, u^ K 0, and c T u = max{0, (1 + 2pp(A)) 1/2 \\ Ab\\ - p P {A)}. 
Proof. Let 

. r i-||A6 /7 ||- 

A = mm|l,pp(A) 

From Proposition 2 it follows that there exists a u such that 

j IIA6 7 | 
Au = A||A6 J ||, u hie 0, \\u\\ < A 



Pp(A) 

Let u = (u,l + Ab HI ,u + Ab 11 , (1 - A)||A6 7 ||,(1 - A)A6 J ). Observe that by con- 
struction 

A T u = b-Ab, u>- K 0. 

Finally, 

c T u= (1 - A)||A& J || 

= max{0, ||A6 7 || + p P (A)(\\Ab n \\ + |A6 m | - 1)} 
< max{0, (1 + 2p P (A))^\\Ab\\ - p P (A)}, 

where the last inequality can be obtained by elementary analysis. □ 

Lemma 11 Let z eAfp. Then if pp(A) > 0, 

a^{^H{w)-^A T ) > (2-28) 

or 

amax(/U 1/2 ^H~ 1/2 ^ T ) < ^ 2 ||#M~ 1/2 II \\A T \\ <2 + 3rp(z). (2.29) 
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Proof. Lemma 9 implies that there exists Ab 6 M m satisfying 



IIA&H < f 



/3 + 2r 



(2.30) 




and such that the optimal value of the following problem 



mm 



c u 



s.t. Au = b + ab, 



(2.31) 



exceeds c^x. Assume that ||A6|| < (1 + 2p 2 p (A))~ 1/2 (c r x + pp{A)). Then by 
Lemma 10 the optimal value of (2.31) is < rf, which contradicts the earlier 
conclusion. Therefore, ||A6|| must satisfy 



3 Finite precision analysis 
3.1 Floating-point arithmetic 

Here we briefly recall the basics of floating-point arithmetic which we will use in 
this paper. A slightly more extensive introduction is in [5, §7]. Detailed treatments 
can be found in books on numerical linear algebra such as [6]. 

We call floating-point numbers a set Fcl containing 0, rounding map a trans- 
formation round : H — > W and round-off unit a constant iieE satisfying < u < 1. 
We require for such a triple that the following properties hold: 

(i) For any x G F, round(x) = x. In particular round(O) = 0. 

(ii) For any i£R, round(x) = x(l + S) with \S\ < u. 

We also define on F arithmetic operations following the scheme 



\\Ab\\>(l + 2p P (A)r^(c T x + PP (A)). 



(2.32) 



Putting (2.30) and (2.32) together, we get 



(p 1/2 H(w)' 1/2 A T ) > 



(l-lf3)(Fx + p P (A)) 
(/? + 2r)(l + 2p 2 (A))V2- 




xoy = round(a; o y) 
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for any x,y G F and o e { + , — , x,/} so that 

o : F x F ->• F. 

It follows from (ii) above that, for any x, y G F we have 

xoy = {x o y)(l + 5), \8\<u. 

We will also use a floating-point version of the square root which, similarly, 
satisfies 

= ^(1 + J), |<5| < u . 

When combining many operations in floating-point arithmetic, quantities such as 
nr=i(l + &i) Pl naturally appear. The proof of the following propositions can be 
found in Chapter 3 of [6]. The notation they introduce, the quantities 7 ra and 6 n , 
and the relations showed therein, will be widely used in our round-off analysis. 

Proposition 3 If \8i\ <u, pi G { — 1, 1} and nu < 1 then 

n 

+ = i + o n 

1=1 

wiiere 

nu 

Wn\ < In 



1 — nu 

□ 

Proposition 4 For any positive integer k such that ku < 1 let Ok be any quantity 
satisfying 

ku 

\0k\ < Ik 



1 — ku 

The following relations hold. 
1) (l + fc )(l + 0j) = l + ^+j, 
2) 

i + Qk _ { i + e k+j ifj<k 

l + 9j \ l + k +2j ifj > k, 

3) If ku,ju < 1/2 then -y k jj < 7 min { feJ }, 

4) ?7fc < 7ifc, 

5) Jk + u< 7fc+i, 

6) 7fc +7j + 7fc7j < 7fc+j- D 
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When computing an arithmetic expression q with a round-off algorithm, errors 
will accumulate and we will obtain another quantity which, we recall, we denote by 
f 1(g). We will also write Error(g) = \q — fl(q)\. 

An example of round-off analysis which will be useful in the sequel is given in 
the next proposition whose proof can be found in Section 3.1 of [6]. 

Proposition 5 There is a round-off algorithm which, with input x,y G IR n , com- 
putes the dot product of x and y. The computed value fl((x,y)) satisfies 

fl((x,y)) = (x,y) +9\i og2 „-]+i(\x\,\y\} 

where \x\ = (\x\\, . . . , \x n \). In particular, if x = y the algorithm computes f l(||x|| 2 ) 
satisfying 

fl(||x|| 2 ) = ||x|| 2 (l + ^ log2 „ 1+1 ). □ 

The following result deals with summation errors. The proof can be found in 
[6], Section 4.2. 

Proposition 6 There is a round-off algorithm which, with input x G H™, computes 
the sum of X{. The computed value f x i) satisfies 

(n \ n n 

i=l J i=l i=l 

In the next section we will have to deal with square roots. The following result 
will help us to do so. 



Proposition 7 Let 6 G R such that \0\ < 1/2. Then, ^1 + 9 = 1+9' with \9'\ < \9\. 
In particular, for a > 

f 1 (Va(l + ^)) = v^(l + 9 k+1 ). (3.33) 

Proof. By the intermediate value theorem we have that y/1 + 9 — 1 = I^Kv^f)' 
with f € (1 - |0|, 1) if < 0, £ € (1, 1 + 9) otherwise. But 



2v^ 



1 



the last since |^| > 1/2. 

Then (3.33) follows from the above. □ 

Our choice of u = 4>(u,(w)), for the function (f> in (3.34) below, guarantees that 
ku < 1/2 holds whenever we encounter 9 k , and consequently, 9 k < 2ku. We will 
therefore not bother the reader by repeating this fact each time we use it. 
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3.2 The finite precision algorithm 

In this section we present a finite precision algorithm that determines which one of 
(P) or (D) is strictly feasible and provides a solution. In the case when the dual 
problem (D) is feasible, after sufficiently refining the precision we will be able to 
obtain an exact feasible solution to (D), however, for the primal problem only an 
approximation to a feasible solution is possible due to the structure of the problem: 
we cannot compute a point on the linear subspace Ax = exactly with finite 
precision. However, we can obtain a forward-approximate primal solution of any 
desired accuracy. To describe this in more detail, we need the following definition 
of a 7-approximate solution. 

Definition 2 Let 7 € (0, 1). A point x G IR n is a ^-forward solution of the system 
Ax = 0, x >-k 0, if x >~k 0, and there exists x 6 M n such that 

Ax = 0, x y K 

and 

\\x — x\\ < tII^H- 

The point x is said to be an associated solution for x. A point is a forward- 
approximate solution of Ax = 0, x >zk 0, if it is a 7- forward solution of the system 
for some 7 € (0, 1). 

We are now ready to present our main result and give a precise description of 
the related algorithm. The proof of Theorem 1 is deferred to Section 4. 

Theorem 1 There exists a finite precision algorithm which, with input a matrix 
A £ jn mxn anc l a number 7 £ (0, 1), finds either a strict ^-forward solution x £ H n 
of Ax = 0, x 0, or a strict solution y £ R m of the system A y <k 0. The 
machine precision varies during the execution of the algorithm. The finest required 
precision is 

u* = (c{n + mfl 2 r 11 *C{A) 7 l 2 )~\ 

where c is a universal constant. The number of main (interior-point) iterations of 
the algorithm is bounded by 

O (rV^logM + \og(C(A)) + I log 7 |)) 

if (P) is strictly feasible and by the same expression without the \ logj\ term if (D) 
is. 

Remark 1 In the numerical analysis literature, fixed precision is used more com- 
monly than variable precision. We note here that from our variable precision analysis 
we can obtain a fixed precision one. Indeed, assume the precision u is fixed. Then 
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our algorithm could run until the point at which it should get a precision finer than 
u. If it found the answer before this point it could return it (and this answer would 
be guaranteed to be correct). If not, it could halt and return a failure message. Fur- 
thermore, the only reason for u to be insufficient is that C(A) is too large. Solving 
the bound for u in Theorem 1 we obtain a lower bound C u for C(A). Thus, the 
failure message could be something like "The condition of the data is larger than 
C u . To solve the problem I need more precision." 

We are now ready to describe our primal-dual algorithm. This is essentially 
a extension of Algorithm IP from Section 2 with some additional features. One 
of these features is the stopping criteria and the other one is the presence of finite 
precision and the adjustment of this precision as the algorithm progresses. To ensure 
the correctness of the algorithm, the precision will be set to 

<KM*)) : = , I ^ Z)1 ' 2 \TT7n (3-34) 

^ VPV JJ crn 5 /2(2 r/ u + l)H/2 v 1 

at each iteration. Here c is a universal constant. 
Let ft = and S = ^. 

Algorithm FP (A, 7) 

(i) Set the machine precision to u := — 77 — — tft, 

\ / 1 cr'(m+n) J /^ 

z := (ae, 1, ae, 2M, -aAe, 0, M e - M m m _m 1 ^ 

\ ' ' ' ' ' ' a ' or ' a ' or ' oc ' ' / 

(ii) Set the machine precision to u := (/)(/j,(z)). 

(iii) If for i = 1, . . . ,r 

Sio - - 6(i(z)r > 

then HALT and return y as a strict solution for A T y -<k 0. 

(iv) If a niill (H(x)-^ 2 A T ) > then HALT and 
return x as a 7-forward solution for Ax = 0, x 0. 

(v) Set]J:= (l-^)MW- 

(vi) Update z by solving the linearization (2.3) of (2.2) for [i = ~fi. 

(vii) Go to (ii). 

The matrix H{w) used in step (iv) is the upper-left nx n block of H(w), where 
w is the scaling point of (x, s) . 

The precise way we solve the system in (vi) is as follows: 

(a) Compute a solution Ay of 

{AH(w)- 1 A T )Ay = AH(w)- 1 (s+JIg(x)). (3.35) 
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(b) Let y:=y + Ay, 



( t ) := ( O" 5 ) m^T^Ay- Ofc(I) + x)). 



Then set x := x + Ax and r := r + Ar. 
(c) Let 



x' := x s' := — y' 

t := 1 t s :=rj 

x" := —Ax s" := — y 
s := 



r, = 1. 



Remark 2 The finite-precision errors in the computations in (b) and (c) are neg- 
ligible compared to the errors involved in solving the linear system on step (a). 
Therefore, for ease of exposition, we will assume that the computations in (b) and 
(c) is step (vi) are exact. We also assume that the initial point z in step (i) and the 
value of ~p in step (v) of Algorithm FP are computed exactly. We stress that these 
assumptions have no consequences in the complexity or accuracy bounds. By mak- 
ing them we can greatly reduce the length of our exposition and focus our analysis 
on the critical stages of the algorithm. 

Under the assumption of infinite precision on steps (b) and (c) the next point 
z + thus defined lies in the linear subspace {Ax = b, A y — s = c}. Moreover, 
Az = (Ax, Ay, Az) satisfies system (2.4) for some (possibly large) r. 

The crux of our finite precision analysis is the estimation of the floating-point 
errors in step (a) above, which we present in Section 3.3. That analysis relies on 
the following technical lemma. 

Lemma 12 Assume z £ Mr and let w = w(z) he its scaling point. With precision 
u = (f)(u,(z)) we can compute B = H(w)~ 1/2 A rr and D = H(w)~ 1/2 A T (where H{w) 
is the upper-left n x n block of H(w) ) satisfying 

1 u 2 

\\i\(B) - B\\,\\fl(D) - D\\ < 7 —!t ^ (3.36) 

11 v ' "' 11 v ' 11 ~ 336-240 r(2r/x + l) 2 v ' 

as well as q = —H(w)~ 1,2 (jlg(x) + s) satisfying 

1 u 3/2 

»"<«> " 0» S 1^240 ' < 3 ' 37 > 
The proof of Lemma 12 in turn relies on the following technical result. 
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Lemma 13 Let z € Mb, and the finite-precision computations are performed with 
u = <f>(n(z)). Then 

( detsj \ detSj 

Error < 7m, Error (detSjdetXj) < det s^det x^m, (3.38) 

\detxij detXi 

where 

(2r^ + l) 2 (log 2 (m + n) + 2) 
^(1-/3)2 

Proof. Observe that from Proposition 5 

Error(detSi) = ||si|| 2 r i og2n .] +1 , Error(det a;,) = Hxill^iog^l+i (3.39) 
for all % G {1, . . . , r}. Let Ki := [log 2 ni\ +1. From (3.39) we have 

det Si 



Error 



/detsj\ det Sj + ||sj|| 2 7 K - , 

< - - ' (1 + ji 

\detx,J detxj - ||xi|| 2 7 Kj 



and from Lemma 7 



det Xj 

_ det Si 7idets^detXi + ||gj|| 2 7 K . + idet Xj + ||xi|| 2 7 Ki det Sj 

detXi (detxj - ||xj|| 2 7 K JdetSi 

< 2 ^tsj f \\xi\\ 2 ||^|| 2 \ 

~~ detxj V detXi Kl det Sj Kl+1 / ' 



* U 112 ' uu °i — II ||2 

W^i II 



Therefore, using Lemma 4 



Error I dets ^ < dets * ll^lPlkipT/ci + lki|| 2 ||*i|| 2 7«i+i + 7i 



detXi J detXi 2// 2 (l - /3) 2 

det Si (2r/x + 1) 2 7ki + (2r/i + 1) 2 7ki +i + 7i 



< 



detXi 2/x 2 (l-/3) 2 



det gj (2r/u + l) 2 
" detx, ' ^ 2 (l-/3) 27Kl+1 ' 

which yields the first inequality in (3.38). The second relation is obtained analo- 
gously. □ 

Proof of Lemma 12. It is well-known (see [17, §3.2]) that the scaling matrix 
H has a block-diagonal structure, where each block corresponds to a Lorentz cone; 
moreover, each individual block can be represented as follows 



Hi{ W { Xi ,Si))- l l 2 = 1 ~ 1 



a -C T 

"» 1 ^ l+a 
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where 7 



*iO II -a. II 



1/4 



a = detfg) 1 / 2 ' and C = det(V/ 2 With C = = 
(7 _ls i0 + T^io, 7~ 1 «7 - 7x7) for alH = 1, . . . , r. 

From Proposition 5 for all i = l,...,r, Error ( (xj, Sj)) < ||xj|| ||si||7i g 2 (n+m)+i> 
and using Lemmas 4 and 7 

Error ((xi, s^) < ( Xi , Si ) p) TW"+™)+i- 

Observe that 

det£ = 2 (\/det x;det Sj + (sj,Xj)j , 

hence by Proposition 7 and Lemma 13 Error(det^) < det^7M+3- 
Further 

Error(7 _1 a) = Error ^ 7 J^_L — ^ = 7~ 1 q6> 3 m+4; (3.40) 

hence 

Error((7 _1 C) i ) < 7~ 1 a7 3 M+4- 

It remains to evaluate the errors in the bottom-left block of Hi(w(xi, Sj)). We have 
Error (^ fc ) < &72M+4, then 

( CfcO \ P f 06 \ . £ 2 Co 2 



„ ( -1 CfcO \ ^ -1 Co 

Error 7 — — < 7 — — 711^+33; 
\ 1 + a J I + a 

-1 (k(l \ ^ __i , Co 



and 



Error + 7 -_J < 7 "^1 + 7 HM + 33. 

Finally, we have 

||f l(H(w)~ 1/2 ) - H( W y 1/2 \\ < (m + lfr" 1 max |a 7 3M+4, (l + ^f^J HlM+u} • 
Observe that by Lemma 7 



det£ = 2(A/detXidetSi + (xi,Si)) > 8(1 - 
Co = 7-^0 + 7x0 < JA= ; 7^< 
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Then 



dete /2 -4^(1-/3)' 1 ] 

Co i i Co ^ 1 i a ^ i i x + s 



i + 7-7^ — = 1 + , ; u t < 1 + Co < 1 + 



l + a detC + Co" 4^(1 - /3) ' 

Observe that Xg + Sg < 2(2r/i + l) 2 . Hence we have 

-if f Co 2 \ 1 5(2r/i + l) 7 /2 
7 max |a73M+4, I 1 + 1 7; ~ J 711M+34 J < 77^ 711M+34- 

Therefore, 

5(n + m)(2r/i + l) 7 / 2 



||fl(tf M~ 1/2 ) - #(w) _1/2 || < ' 711M+34 

and 



Ai 3 /2 



llfl^H" 1 ^) -^H- 1/2 || < 1)7/2 711M+34- 

Now we estimate the error in computing B = H(w)~ 1 ^ 2 A T . First, observe that the 
error for multiplication of f l(H(w)~ 1/2 ) by A T can be estimated as follows (see [6, 
Chapter 22]) 

||fl(iJH- 1/2 )^ T - fl (fl(H(w)- 1/2 )A T ) || < n«||fl(fr(«;)- 1 / 2 )||||^ T ||. 
Therefore 

\\B - fl(B)\\ < \\H(w)- 1/2 A T - fl(iT(w)- 1/2 ).A T || 

+ \\fl(H(w)- 1/2 )A T - f 1 (fl(H(w)~ 1/2 )A T ) || 

< (\\H(w)- 1/2 - fl(H( W y i/2 )\\ + nu\\fl(H(w)- 1/2 )\\) \\A T \\ 
20(n + m)(2r/u + I) 7 / 2 

< -772 712M+34- 

Since our precision u satisfies (3.34), this yields (3.36) for B. The corresponding 
bound for D is obtained analogously. 

It remains to evaluate the errors in computing q. By straightforward computa- 
tion we obtain for each 'block' 

q t0 = -het£ /2 (l-n 1 ' 2 



det Xi 



-2 

i 

det Xi 



Ql= "dete 1/2 + 6 



Sio + -7idetC- /2 ) Xi + ( -7^ ^et^ 2 + x i0 I s; 



HO 

We obtain (3.37) by a similar argument as when evaluating ||f l(H(w)~ 1/2 ) — 
H(w)~ 1/2 \\. We omit this tedious exercise for the sake of brevity. □ 
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3.3 Finite-precision analysis of solving the Newton system 



The main result of this section is Lemma 14, which bounds the round-off error in 
the computation in the reduced equations (3.35) in (a), when it is performed with 
finite precision. 

Observe that the reduced system of equations (3.35) is equivalent to the least- 
squares problem 

min \\Bv + q\\ 2 

for B = H(w)~ 1 / 2 A T and q = —H{w)~ 1/2 (jlg(x) + s). We rely on this equivalence 
to obtain the bound in Lemma 14. More precisely, we apply a known round-off 
error result for the least-squares problem, namely Proposition 8. Lemma 14 follows 
from Proposition 8 and suitable bounds on the norm of q and singular values of B 
obtained earlier in Lemmas 6 and 11 respectively. 

Recall the following stability property of Golub's method for linear least-squares 
(cf. [7, Chapter 16]). 

Proposition 8 Let B € W xl (p > I) have full rank. Let u denote the machine 
precision. If Golub's method is applied to 

min \\Bv + /|| 

vGTRl 

the computed solution is the exact solution to a problem 

mm \\(B + 8B)v + (f + 6 f)\\ 

where 

\\SB\\ < cupl 3 / 2 \\B\l 11,5/11 < cupl || /|| 
and c is a universal constant independent of p and I . □ 



Lemma 14 Let z G Np. With precision u = 0(/x) in all arithmetic operations, we 
can compute a vector f l(Ay) such that 



\\{AH{w)-'A L )±l{&y)- AH{w)-\plg{x) + s))\\ < 



240r(2r^(z) + 1)' 
In addition, fl(Ay) < 12r/i~ 1/2 . 

Proof. Let B = H(w)~ 1/2 A T and q = H(w)~ 1/2 (j[g(x) + s), then the system of 
equations 

(AH{w)- 1 A T )Ay = AH{w)- 1 (jig(x) + s) 

can be written as B T BAy = —B T q and its solutions are those of the least squares 
problem 

min ||Bv + g||. (3.43) 
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Hence we can apply Golub's method to (3.43) to compute a solution to the original 
equation. Let Ay be the vector actually computed by Golub's method when solving 
(3.43). Then Ay is the exact solution of 

min \\Bv + q\\ (3.44) 

ti£E m 

for some B and q satisfying 

\\B - f 1(S)|| < cum 3/2 n||f l(B)||, \\q- fl(g)|| < cumn||f l(g)||, (3.45) 

where c is a universal constant. Let AB = B — B, Aq = q — q. Since Ay is an exact 
solution of the least squares problem (3.44), we have B T BAy + B T q = 0, and thus 

B T BAy + B T q = -(B T AB + AB T B)Ay — (B T Aq + AB T q). (3.46) 

From Lemmas 5 and 12 and ||^4 T || < y/2 we have 

||fl(B)|| < || J ffH- 1/2 ||P T || + ||fl(5) - B\\ < 3(2l "^ + 1) . (3.47) 
Analogously, Lemmas 6 and 12 yield 

iifi(,)ii<Ni + iifi(,)- 9 ii< V + (2r7Tiy- (3 - 48) 

Then from (3.45), (3.47), (3.48) and our choice of u we have 

\\B - f 1(B) || < 3cum 3 / 2 n^±I < 1 fl ■ (3.49) 

11 V 711 - ^ 1/2 - 336 . 480 r (2 r/ u + l) 2 ' v ' 

1 u 3 / 2 

\\q- i\(q)\\ < 3cumn/i 1/2 < — ^ -. (3.50) 

" y W " ~ P " 16-480 r(2r/i + l) 2 V ; 

Applying Lemma 12 again and using (3.49) and (3.50), 

1 u 2 

\\AB\\ < \\fl(B) - B\\ + \\B - f 1(5)11 < -7- ; (3.51) 

ll ll - ll v. ) 11 11 v ^"- 168 .48o r (2r^ + l) 3 / 2 ' V ; 

1 u 3 / 2 

l|A«|| < IHK,) - „|| + |f - f 1(„)|| < 5^5 ■ (3.52) 
Then from (3.51), (3.52) and using the bounds on ||B|| and \\q\\ discussed above, 

||B|| < + \\AB\\ < 4^^; \\q\\ < \\q\\ + \\Aq\\ < ^\ (3.53) 

From (3.51) and (3.53) we have 



1 11 3 / 2 

\B T AB + AB T B\\ < (\\B\\ + ||AB||)||AB|| < f -. (3.54) 

1 t 11 - vii 11 -r 11 11/11 H- 480 I2r(2r/i + l) V 1 
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From (3.51), (3.52) and (3.53) 

\\B T Aq + AB T q\\ < ||£||||Ag|| + ||AB||||g|| < ^ • (3.55) 

It remains to bound ||Ay||. Using Lemma 11 and (3.51) we have 
a min (B) > a min (B) - \\AB\\ > HL- JL = JL. 

minV ) - mini ) || || _ ^ ^ ^ 

From (3.53) we have \\q\\ < fi 1/2 . Observe that since \\BAy + q\\ = min„ \\Bv + q\\, 

\\Ay\\ < —S-^ < Yl^l\ (3.56) 

Finally, we have from (3.46), (3.54), (3.55) and (3.56) 

\B T BAy + B T q\\ < \\B T AB + AS T B|| • || Ay|| + \\B T Aq + AB T q\\ 
1 M 3/2 1/2 1 fi 



480 12r(2r/i + l) ^ 480 2r/x + 1 

240r(2r/i + 1) 



< -A r. □ 



3.4 Finite-precision analysis of termination conditions 

Lemma 15 (Dual termination) Let Pd(A) > and z G Np with fi(z) < ■ 
Then 

f 1 (s i0 - \\sl\\) > f 1 (Gfi(z)r) for i € 1 : r. (3.57) 

Moreover, if z G Np satisfies (3.57), then the subcomponent y of z is a strict feasible 
solution to (D); in other words, A T y -<k 0. 

PROOF. From our choice of precision u = (j)(fi(z)) and the fact that z G Afp it 
readily follows that 

Error(s i0 - < Tfi(z) (3.58) 

On the other hand, by Lemma 3 we have 

Sl ° ~ ^ ~ \r~^ PDiA) ~ 27v^ (4 ° r3/U(z)) " 10 ^ r - 

Thus 

f l(s i0 - \\si\\) > f l(8r/i(z)) > f l(6r//(z)). 
Now assume (3.57) holds. Again, by (3.58) we get 

Sio ~ \\s~i\\ > 5rfj,(z). 
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Since Ay + s = c, in particular A T y — y' + s = 0. So —A T y = s — y'. Since 
||y'|| < ?7 < r + V = c^x — t^y = 2r/x(z), it follows that for i = 1 : r we have 

s *o - y'io ~ \\s~i ~ yiW > s i0 - \\sl\\ - y' i0 - \\y-W > 5rn.(z) - 2\\y'\\ > 17/(2:) > 0. 



Lemma 16 (Primal termination) Assume z G Mp- If /j.(z) < + ^ 

then in step (iv) the algorithm yields 



Moreover, if z £ Mp satishes (3.59) then the subcomponent x of z is a ^-forward 
solution of Ax = 0, x 0, and 



is an associated solution for x. 

Proof. Let D = H(x)~ 1/2 A T and assume that we compute <J m i n (D) using a 
backward stable algorithm (e.g., QR factorization). Then the computed f l(<j m i n (D)) 
is the exact <7 m i n (fl(D) + E) for a matrix E with ||D|| < cn 2 u\\i 1(D) || for some 
universal constant c (see [6]). We have 



Therefore, s — y' >~k and consequently A y = y' — s <k 0. 



□ 





(3.59) 



x = x- H(x)- 1 A T (AH(x)' 1 A T )- 1 Ax 



Error(<7 min (D)) < ||f 1(D) - D|| + cn 2 u(\\D\\ + ||f 1(D) - D||). 



(3.60) 



By our choice of u and bounds from Lemmas 5 and 12, we get 



Error(cr min (D)) < 17/(2). 



Therefore 



f l(a min (D)) > 

^min 

(D) - Tfl(z). 



Using the bound from Lemma 2 we have 




By our choice of fx(z) this yields 




Now assume (3.59) holds. Again by (3.60) we get 



<WD) > f 1( "min 



2r [i(z) 



7 
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Denote Ax = -H(xy 1 A T (AH (x) 1 A T ) 1 Ax. From Lemma 4 and \\Ax\\ = \\x"\\ 
we have 

||A,||1 _ (A x f { AH( xr ^r^ < Jgt_ < _gg2L < 7 , 
Furthermore, by Proposition l(a,c) we have 

lj w - wi^i - mx) "^ A = l|Ax|U - 1<l - 

Therefore, x = x + Ax is a 7-forward solution of Ax = 0, x 0. □ 



4 Proof of the main result 

We are finally in a position to prove our main result (Theorem 1). We first prove 
that on every step the algorithm keeps up with the central path, and at the same 
time the value of ~p decreases by a fixed factor. Then we show that once fi is small 
enough to satisfy either dual or primal termination conditions (Steps (iii) and (iv) 
of the algorithm), the algorithm terminates and yields a correct answer. Then the 
bound on the number of iterations follows trivially from the termination bounds on 
[i and the factor of the decrease of [i. Similarly we obtain bounds for the finest 
precision based on the precision update function cj) and the aforementioned bounds 
on fi. 

Before we go ahead with the proof, we need to guarantee that the initial point 
satisfies all the necessary bounds. The following result is an immediate consequence 
of [19, Proposition 4.6]. 

Lemma 17 (Computation of the initial point) The initial point 

( M M M M M 

z := ae, 1, ae, 2M, —aAe, 0, — e, — e, — e, 1, 

\ a a z a a z a 



where a = and M = Q ^ 6 ^ , satisfies z £ Np and fi(z) = a ^ e ^ = 0(1). 

□ 



Proof of Theorem 1. We first disregard the halting steps (iii) and (iv) of the 
algorithm and prove that, no matter how many iterations we have performed, all 
our iterates stay close enough to the central path. We use an induction argument. 

The induction base is given by Lemma 17. We now assume that at the start of 
Step (ii) of the algorithm the value of z satisfies z £ Np. We need to show that the 
vector z + computed in step (vi) is also in Np ■ 
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Prom Lemma 14 we know that the point z + computed with finite precision in 
(a), and infinite precision in (b) and (c), satisfies 

AAx =0 
^l T Ay + As =0 
Ax + H(w)~ 1 As =-(]Ig(s) + x) + r, 

for some r with ||r|| < 12 o r fi+^(z)) • Hence from Lemma 1 we have z + € Np. Hence, 
z € Np on every iteration. 

Now we show the bounds for the number of iterations and the finest precision. 

From Lemma 16 we know that once n(z) reaches the lower bound of 



Pp(A) 



problem. 

It follows from Lemma 1 that 



;0 , the algorithm yields a correct 7- approximate solution to the primal 



120r 2 V y/2r 120r 2 

Therefore, taking into account that r 2 > y/2r ■ 3^/| > 3\/2r and that S = 1/45, we 
have 

H(z + ) < (l ±=) n{z). 



Given an initial value of fi(zo), after k iterations we have 

ti'k) < (l " M(^d). 



Since we want /x(^) < ^1 + ^ , we have the condition 



1 \ k , , WA) / l x " X 

1 ^ )<^y 1+ 



60\/2r/ !0r 2 V 7 

By taking logarithms on both sides, and using fi(zo) = 0(1), we get the desired 
relation 

k = O (rV 2 Qog(r) + log(C(A)) + | log 7 |)) . 

Since the algorithm halts once /i(z) < p [qJ? ^1 + , we deduce from our 
precision update formula (3.34) that the finest required precision u* satisfies 



u* > [c(n + m 

f/\S c{A) 7/2\ 
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Similarly, for the dual feasible case from Lemma 15 we know that /j,(z) < p ^ r3 
guarantees successful termination of the algorithm. Hence we similarly get the 
bound 

k = O (V/ 2 (log(r) + log(C(A)))) , 
and for the finest precision we have 

u* > (c(m + n) 5/2 r 1L5 C(^) 7/2 ) ~* . □ 
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